统计问题｜绘制任意分布的 QQ 图

点击下方公众号，回复资料分享，收获惊喜

简介

最近科研论文的审稿意见，需要对数据的拟合情况进行说明。分位数-分位数图^[1]（Quantile-Quantile plot，简称 QQ 图）是一种不错的选择。

研究动机

在 R 语言中常用 qqnorm() 绘制 QQ 图的散点图，qqline() 用于在图上添加一条参考直线。可支持与一些常见的理论分布进行比较，例如：正态分布 ("norm")、指数分布 ("exp")、伽马分布 ("gamma")、威布尔分布 ("weibull")、对数正态分布 ("lnorm")等。

但是绘制出来的图形可能不是很美观，如下所示：

本文贡献点

当然小编以前使用 ggplot 绘制过：绘制不同分布的 QQ 图，但内容不包括逆高斯分布以及其他比较“新”的分布。

基于这个痛点，本文小编将根据 QQ 图的定义，通过 ggplot2^[2] 复现自己常用的三个分布（正态、伽马和逆高斯分布）的 QQ 图。

题外话：读者可以根据本文思路，拓展到任何自己感兴趣的分布。

教程

理论介绍

介绍：一种用于检查数据是否符合某个理论分布的统计图。

基本思想：比较观察数据的分位数与理论分布的分位数。

判断依据：如果数据符合理论分布，那么 QQ 图上的点将近似落在一条直线上。

绘制步骤：

1. 排序数据：将观察数据按照大小进行排序。

2. 计算分位数：计算每个观察数据点对应的分位数。这可以通过排名除以总数据点的数量来实现。

3. 计算理论分位数：使用所选择的理论分布的累积分布函数，计算对应于每个观察数据点的理论分位数。

4. 绘制 QQ 图：将观察数据的分位数作为横坐标，理论分位数作为纵坐标，绘制散点图。如果数据符合理论分布，散点图将近似为一条直线。

代码实现

根据上面的理论结果，我们依次进行代码实现。假设数据来自于 gamma 分布。

y <- rgamma(200,10,1)

1. 排序数据

sorted_data <- sort(y)

1. 计算样本分位数

sample_quantiles <- (1:length(sorted_data) - 0.5) / length(sorted_data)

1. 计算理论分位数

首先，对数据进行参数估计，得到估计之后，产生理论分位数。

3.1 参数估计

这里使用极大似然方法估计参数，首先定义对数似然函数：

library(statmod)
fn_ig <- function(par= c(mean, dispersion), data) {
  ll <- -sum(dinvgauss(data, mean= par[1], dispersion= par[2], log= TRUE))   
  return(ll)
}

3.2 参数估计，并计算理论分位数：

pp <- optim(par= c(1,1), fn_ig, data = y)
mq <- pp$par[1]; disp <- pp$par[2] 
theoretical_quantiles <- qinvgauss(sample_quantiles, mean= mq, dispersion= disp)

1. 绘制 QQ 图

数据整理并绘制

#数据整理
data_ig = data.frame("x" = theoretical_quantiles, "y" = sorted_data)
#绘制
library(ggplot2)
ggplot(data_ig,aes(x=x, y=y)) +
  geom_point(color="#4CA6A3",size=1.2) + 
  geom_abline(intercept = 0, slope = 1,color="#693476",size=1) + 
  # facet_wrap(vars(factor(pc)),scales = "free") + 
  theme_bw() + theme(legend.position = "none") + 
  xlab("Theoretical") + ylab("Sample")

此时，给出了该数据使用逆高斯分布拟合的情况，结果还算可以！读者可以通过一些检验（Kolmogorov–Smirno^[3]等），进行进一步检测。

拓展

根据这个思路，小编写了正态和 Gamma 的结果，并通过一个函数将三种分布的整合在一起。最后绘制出三种分布拟合该数据的结果。

还计算了三种分布拟合的 KS 检验的 p 值：

p 值均大于 0.05，接受原假设，即三种分布拟合该数据均可。读者可以将合并的图放到自己论文中，仅需将 y 修改成自己的数据即可！

汇总代码

y <- rgamma(200,10,1)
library(statmod)
fn_ig <- function(par= c(mean, dispersion), data) {
  ll <- -sum(dinvgauss(data, mean= par[1], dispersion= par[2], log= TRUE))   
  return(ll)
}
fn_normal <- function(par = c(mean, sd), data) {
  # 正态分布的对数似然函数
  ll <- -sum(dnorm(data, mean = par[1], sd = par[2], log = TRUE))
  return(ll)
}
fn_gamma <- function(par = c(shape, rate), data) {
  # 伽马分布的对数似然函数
  ll <- -sum(dgamma(data, shape = par[1], rate = par[2], log = TRUE))
  return(ll)
}
qqdat = function(den = "ig", y = diff_y_all, par= c(1, 1)){
  # 排序数据
  sorted_data <- sort(y)
  # 计算样本分位数
  sample_quantiles <- (1:length(sorted_data) - 0.5) / length(sorted_data)
  
  if(den == "ig"){
    # 参数估计
    pp <- optim(par= par, fn_ig, data = y)
    mq <- pp$par[1]; disp <- pp$par[2] 
    # 计算理论分位数
    theoretical_quantiles <- qinvgauss(sample_quantiles, mean= mq, dispersion= disp)
  } else if(den == "norm"){
    pp <- optim(par= par, fn_normal, data = y)
    mq <- pp$par[1]; disp <- pp$par[2] 
    theoretical_quantiles <- qnorm(sample_quantiles, mq, disp)
  } else if(den == "gamma"){
    pp <- optim(par= par, fn_gamma, data = y)
    mq <- pp$par[1]; disp <- pp$par[2] 
    theoretical_quantiles <- qgamma(sample_quantiles, shape = mq, rate = disp)
  }
  
  data_ig = data.frame("x" = theoretical_quantiles, "y" = sorted_data)
  return(list(data_ig,"Para"= c(mq,disp)))
}

# 数据处理
qqplot_den = function(den = "norm",data6 = y,par= c(10, 1)){
  u1 = qqdat(den = den, data6, par= par)
  ig_g = data.frame(u1[[1]])
  
  w1 = ggplot(ig_g,aes(x=x, y=y)) +
    geom_point(color="#4CA6A3",size=1.2) + 
    geom_abline(intercept = 0, slope = 1,color="#693476",size=1) + 
    theme_bw() + theme(legend.position = "none") + 
    xlab("Theoretical") + ylab("Sample")
  
  if(den=="norm"){
    kstest = ks.test(data6, "pnorm", mean = u1$Para[1], u1$Para[2])$p.value
  } else if(den=="ig"){
    library(fitdistrplus)
    kstest = ks.test(data6, "pinvgauss", mean = u1$Para[1], dispersion = u1$Para[2])$p.value
  } else if(den=="gamma"){
    kstest = ks.test(data6, "pgamma", u1$Para[1], u1$Para[2])$p.value
  }
  
  return(list(w1,kstest))
}

w1 = qqplot_den(den = "norm",data6 = y,par= c(10, 1))
w2 = qqplot_den(den = "ig",data6 =  y,par= c(1, 1))
w3 = qqplot_den(den = "gamma",data6 =  y,par= c(1, 1))
cowplot::plot_grid(w1[[1]],w2[[1]],w3[[1]],labels=c("N","IG","G"),nrow=1)
# ggsave("qqplot.pdf", width = 8,height = 4)

参考资料

[1]

分位数-分位数图: https://haomenghit.github.io/2019/08/04/Q-Q%E5%9B%BE%E4%BB%8B%E7%BB%8D/

[2]

ggplot2: https://ggplot2.tidyverse.org/

[3]

Kolmogorov–Smirno: https://en.wikipedia.org/wiki/Kolmogorov%E2%80%93Smirnov_test

相关推荐